% find out the unfolding steps automatically for constant loading rate
% experiment. Unfolding force is from magnet position; while step size is
% calculated based on two linear fitting line before and after the
% unfolding transition (help: Linear Least Squares).

% to do: remove from the one with smallest t-value and one by one.
% key method: increase t_threshold gradually, then no need to remove from smallest t-value.

clc; 
clear all;
close all;

%%  path and file name w/o '.txt'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20110808-Cu-Na150-M270\B4'  %%%%
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20110808-Cu-Na150-M270\B8'  %%%%
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20110812-Cu-Na150\ch2-B1'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20110812-Cu-Na150\ch2-B5-BSA'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20110831\B1'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20110831\B6'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20110831\ch2-B1'
%  name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20110831\ch2-B13'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20110831\ch2-B14'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20110831\ch2-B4'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20110831\ch2-B8'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20111123-Cu-Na150-pure-bf\Ch2-B1'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120413-2.1pNpsec-Cu-BSA-Na150\B4'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120413-2.1pNpsec-Cu-BSA-Na150\B5-morebd'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120413-2.1pNpsec-Cu-BSA-Na150\B6-morebd'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120413-2.1pNpsec-Cu-BSA-Na150\ch2-B10-morebd-3'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120413-2.1pNpsec-Cu-BSA-Na150\ch2-B9-morebd-3'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120413-2.1pNpsec-Cu-BSA-Na150\ch3-B1'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120414-2.1pNpsec-Cu-BSA-Na150\B14-moreFLNA-morebd3'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120414-2.1pNpsec-Cu-BSA-Na150\B2'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120414-2.1pNpsec-Cu-BSA-Na150\B3'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120414-2.1pNpsec-Cu-BSA-Na150\B5-morebd'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120414-2.1pNpsec-Cu-BSA-Na150\B9-morebd2'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120414-2.1pNpsec-Cu-BSA-Na150\ch2-B3'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120521-Cu-Na150\B1-BSA'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120521-Cu-Na150\B7-BSA-morebd'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\wang yejun\IgFLNa1-3-20-21-1-3\20110909\channel2\B4'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\wang yejun\IgFLNa1-3-20-21-1-3\20110909\channel3\B3'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\XiaoMeng\2012-04-11\20-21_Ch1_B14'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\XiaoMeng\2012-04-11\20-21_Ch2_B2'

%% 0.1 pN/s
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120323-Cu-BSA-Na150-M270-LED\B4-morebd'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\IgFLNa1-3-20-21-1-3\20120324-Cu-BSA-Na150-M270-LED\B5-0.1pNpsec'
 name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\Saranya\IgFLNa20-21\20120321-BSA-Na 150-M270-pump\B4-morebd-0.1pNpsec'
% name='C:\Users\mbich\ChenHu74\RESEARCH\Experiment\IgFLNa20-21\Saranya\IgFLNa20-21\20120321-BSA-Na 150-M270-pump\ch2-B5-0.1pNpsec-morebd'
% 
% ---- find the maximal d_ext values and locations ----
% magnet > 11, d_ext > a certain threshold, then find the local maximal of
% d_ext, that's the transition point.
magnet_threshold = 10; %9;  %%%%
step_threshold = 3; %10; %%%%
t_threshold_final = 10;  %%%% %%%%  9
num_std = 4; %2; %3; %4; %%%%%
lifetime_threshold = 2;

smooth_span = 10; %%%% Size of the averaging window

%% input data file
filename=sprintf('%s%s',name,'.txt');
rawdata=importdata(filename,'\t',2);
mydata=rawdata.data;

%get individual columns
time=mydata(:,1);
magnet=mydata(:,2);
ext=mydata(:,3);
dx=mydata(:,4);
dy=mydata(:,5);
N=size(time,1);

% smooth the extension
smooth_ext = medfilt1(ext,smooth_span);  % median keeps the sharp transition better, but will remove small peaks
% smooth_ext= smooth(ext,smooth_span,'sgolay',2);  % also good to keep sharp transition.
subplot(1,2,1);
plot(time,ext,'b-',time,smooth_ext,'r-');
legend('ext.','Smoothed ext.')
title('Extension time course');
xlabel('Time(sec)');
ylabel('Extension(nm)');

%calculate extension changes locally: d_ext(1:N-span)
span = 2; %4; %%%%  to detect very sharp transtion only, even number
d_ext = smooth_ext(span+1:N) - smooth_ext(1:N-span);
 subplot(1,2,2);
 plot(time(1+span/2:N-span/2),d_ext,'b-')
% plot(magnet(span/2+1:N-span/2),d_ext,'b-')
legend('delta ext.')
title('delta Extension time course');
xlabel('Time(sec)');
ylabel('delta Extension(nm)');

% ---- analyze local fluctuation of d_ext ----
% local variance of d_ext by using guassian distribution fit.
window_var = 100; %%%% one second window to calculate variance.
std_d_ext = zeros(N-span-window_var,1);
ave_d_ext = zeros(N-span-window_var,1);
for i=1:N-span-window_var
    temp = d_ext(i:i+window_var-1);
    [ave_d_ext(i),std_d_ext(i)] = normfit(temp);
    %temp1 = mle(temp);   ave_d_ext(i)=temp1(1); std_d_ext(i)=temp1(2);
end
hold on;
plot(time(1+span/2+window_var/2:N-span/2-window_var/2),ave_d_ext+std_d_ext*3,'g-',time(1+span/2+window_var/2:N-span/2-window_var/2),ave_d_ext-std_d_ext*3,'g-');
hold off;

% exclude peaks when calculating local variance
factor_exclude = 2.3; %%%%  ave+std*factor_exclude, 2.3 for window_var=100, sqrt(log(100))=log(10)=2.3
for j=1:2 %%%% 6, peiwen comments that no more than 1 cycle should be used to exclude outsiders.
    censoring = zeros(N-span,1);
    censoring(1+window_var/2:N-span-window_var/2) = (d_ext(1+window_var/2:N-span-window_var/2)>ave_d_ext+std_d_ext*factor_exclude | d_ext(1+window_var/2:N-span-window_var/2)<ave_d_ext-std_d_ext*factor_exclude);
    keep = ~censoring;  % NOT
    for i=1:N-span-window_var
        temp = d_ext(i:i+window_var-1);
        ave_d_ext(i) = sum(temp.*keep(i:i+window_var-1))/sum(keep(i:i+window_var-1));
        std_d_ext(i) = sqrt(sum(((temp-ave_d_ext(i)).^2).*keep(i:i+window_var-1))/(sum(keep(i:i+window_var-1))-1));
        % [ave_d_ext(i),std_d_ext(i)] = normfit(temp);  % not accurate since the zero elements are also counted.
    end
end
hold on;

plot(time(1+span/2+window_var/2:N-span/2-window_var/2),ave_d_ext+std_d_ext*num_std,'r-',time(1+span/2+window_var/2:N-span/2-window_var/2),ave_d_ext-std_d_ext*num_std,'r-');
hold off;
up_limit = ave_d_ext+std_d_ext*num_std;  % size N-span-window_var
low_limit = ave_d_ext-std_d_ext*num_std;  %

%% trial: find the peak of d_ext(1:N-span)

num = 0;
for i=1:N-span-window_var  % i is index of up_limit, low_limit
    i_ext = i+span/2+window_var/2;
    i_d_ext = i+window_var/2;
    if((magnet(i_ext)>magnet_threshold) && (d_ext(i_d_ext)>up_limit(i) || d_ext(i_d_ext)<low_limit(i)))
        if (d_ext(i_d_ext) == max(d_ext(i_d_ext-span:i_d_ext+span)) || d_ext(i_d_ext) == min(d_ext(i_d_ext-span:i_d_ext+span)) )  % local maximal or minimal in a range.
            if num==0  % the first transtion
                num = num +1;
                unfold_index(num) = i_ext;
                unfold_d_ext(num) = d_ext(i_d_ext);
            elseif i_ext-unfold_index(num) > lifetime_threshold
                num = num +1;
                unfold_index(num) = i_ext;  
                unfold_d_ext(num) = d_ext(i_d_ext);
            else      % exclude short lifetime state
                if abs(d_ext(i_d_ext)) > abs(unfold_d_ext(num))  %replace the old transition with the new one.
                    unfold_index(num) = i_ext;
                    unfold_d_ext(num) = d_ext(i_d_ext);
                end
            end
        end
    end
end
if i==N-span-window_var
    index_end = i_ext-1;
end

%% linear fitting of ext, for linear fitting: y=b1*x+b2, to calculate the step size.
% http://en.wikipedia.org/wiki/Student's_t-test#Unequal_sample_sizes.2C_unequal_variance

span_fitting = 100; %%%% maximal windown size to do linear fitting or average of ext, for linear fitting: y=b1*x+b2
min_step = 0;
min_t = 0;
min_n_right = 0;

for t_threshold=0.1:0.1:t_threshold_final
disp('t_threshold=');
disp(t_threshold);
while min_t < t_threshold || min_step < step_threshold || min_n_right < lifetime_threshold
    step_size = zeros(num,1);
    t_value = zeros(num,1);
    n_left = zeros(num,1);
    n_right = zeros(num,1);
    ext_left = zeros(num,1);
    var_left = zeros(num,1);
    ext_right = zeros(num,1);
    var_right = zeros(num,1);
    for i=1:num
        if i==1
            left = max(unfold_index(i)-span/2-span_fitting,1);
        else
            left = max(unfold_index(i)-span/2-span_fitting, unfold_index(i-1)+1);
        end
        right = unfold_index(i)-span/2;

        x = time(left:right);
        y = ext(left:right);
        n = right-left+1;
        n_left(i) = n;
        if n<=20
            ext_left(i) = mean(y);
            var_left(i) = var(y);
        else
     %       b(1)=(n*sum(x.*y)-sum(x)*sum(y))/(n*sum(x.^2)-(sum(x)^2));
      %      b(2)=(sum(y)-b1*sum(x))/n;
            b=polyfit(x,y,1);
            ext_left(i) = b(1)*time(unfold_index(i))+b(2);
            var_left(i) = var(y-(b(1)*x+b(2)));
        end
     %--------------------------------------------------------------------
        left = unfold_index(i) +span/2;
        if i==num
            right = min(unfold_index(i)+span/2+span_fitting, index_end-1);
        else
            right = min(unfold_index(i)+span/2+span_fitting, unfold_index(i+1)-span/2);
        end    

        x = time(left:right);
        y = ext(left:right);
        n = right-left+1;
        n_right(i) = n;
        if n<=20
            ext_right(i) = mean(y);
            var_right(i) = var(y);
        else
            b=polyfit(x,y,1);
            ext_right(i) = b(1)*time(unfold_index(i))+b(2);
            var_right(i) = var(y-(b(1)*x+b(2)));
        end
    %---------------------------------------------------------------------
        step_size(i) = ext_right(i) - ext_left(i);
        t_value(i) = abs(step_size(i))/sqrt(var_left(i)/n_left(i)+var_right(i)/n_right(i));  % did not consider the correlation length of ext.
    end
    min_step = min(abs(step_size(1:num)));
    min_t = min(t_value(1:num));
    min_n_right = min(n_right(1:num));
    
     time(unfold_index(1:num))
     magnet(unfold_index(1:num))
     step_size(1:num) = step_size(1:num)
     t_value(1:num) = t_value(1:num)
     n_right = n_right
%     figure;
%     subplot(3,1,1);
%     hist(step_size,100);
%     subplot(3,1,2);
%     hist(t_value,100);
%     subplot(3,1,3);
%     hist(n_right,100);
    
    if min_t < t_threshold || min_step < step_threshold || min_n_right < lifetime_threshold  % remove small steps
        num_new = 0;
        for i=1:num
            if t_value(i) > t_threshold && abs(step_size(i)) > step_threshold && n_right(i) > lifetime_threshold  % keep good step
                num_new = num_new +1;
                unfold_index_new(num_new) = unfold_index(i);
%             else  % sometimes necessary, otherwise the num_new can be zero, then it will cause error.
%                 if i<num && n_right(i) < span_fitting  % not the last one, and removing this small step may affect the following step
%                     if t_value(i) > t_value(i+1) || abs(step_size(i)) > abs(step_size(i+1))
%                         num_new = num_new +1;
%                         unfold_index_new(num_new) = unfold_index(i);
%                     end                    
%                 end              
            end
        end
        unfold_index(1:num_new) = unfold_index_new(1:num_new);
        num=num_new;
    end
end
end

%%
filename = sprintf('%s%s',name,'-f-step.dat');
fid = fopen(filename,'w');
if(fid == 0)
    disp('file cannot be opened');
end
fprintf(fid,'num\ttime(sec)\tmagnet(mm)\tstepsize(nm)\text1\text2\tnleft\tnright\tt-value\n');
for i=1:num
    n = unfold_index(i);
    fprintf(fid,'%d\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%d\t%d\t%.1f\n',n,time(n),magnet(n),step_size(i),ext_left(i),ext_right(i),n_left(i),n_right(i),t_value(i));
end
fclose(fid);